Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


A Review of vFDTD: An Efficient and Accurate 
Method to Solve Low Frequency, Non-Conformal 
and Multi-Scale Problems Using FDTD 


Kadappan Panayappan! and Raj Mittra»? 
1 TE Connectivity, USA. ? Electromagnetic Communication Lab, Pennsylvania State 


University, USA. è University of Central Florida, USA. 


rajmittraQ@ieee.org 


Abstract 


In this write up we present a review of vFDTD, an efficient and accurate method to solve 
low frequency, non-conformal and multi-scale problems using FDTD. The conventional 
time domain technique FDTD demands extensive computational resources when solving 
low frequency problems, or when dealing with dispersive media, or when structures with 
high Q factor are involved. The vFDTD (New FDTD) technique is a new general-purpose 
field solver, which is designed to tackle the above issues using some novel approaches, 
which deviate significantly from the legacy methods that only rely on minor modifications 
of the FDTD update algorithm. The yFDTD solver is a hybridized version of the confor- 
mal FDTD (CFDTD), and a novel frequency domain technique called the Dipole Moment 
Approach (DM Approach). This blend of time domain and frequency domain techniques 
empowers the solver with the potential to solve problems that involve: (i) calculating 
low-frequency response accurately and numerically efficiently; (ii) handling non-Cartesian 
geometries such as curved surfaces accurately without staircasing; (iii) handling thin struc- 
tures, with or without finite losses; and (iv) dealing with multi-scale geometries. 
Keywords: Method of Moments (MoM), Finite Element Method (FEM), 
Finite Difference Time Domain (FDTD) Method, Dipole Moment Approach 
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(DM Approach), Low Frequency (LF) 


0.1 Introduction 


The conventional time domain technique FDTD (Finite Difference in the Time Domain) 
demands extensive computational resources when solving low frequency problems, or when 
dealing with dispersive media or when structures with high Q factor are involved. To 
tackle some of these challenges, the conventional techniques are often modified in a manner 
that is tailored to solve a particular problem of interest. However, more often than not, 
these tailored methods turn out to be computationally expensive, and they often lead 
to instabilities. Hence, it is useful to develop techniques that can overcome the above 
limitations, while preserving the advantages of the existing methods. The vFDTD (New 
FDTD) technique, which is described in this chapter, is a new general-purpose field solver, 
which is designed to tackle the above issues using some novel approaches, that deviate 
significantly from the legacy methods that only rely on minor modifications of the FDTD 


update algorithm. 


0.2 vEFDTD Solver 


The vFDTD solver is a hybridized version of conformal FDTD (CFDTD) [1], and a novel 
frequency domain technique called the Dipole Moment Approach (DM Approach) described 
in [2, 3]. This blend of time domain and frequency domain techniques empowers the solver 


with the potential to solve problems that require: 
e Calculating low frequency response accurately and numerically efficiently 


e Handling non-Cartesian geometries such as curved surfaces (see Fig. 1) accurately 


without staircasing 


e Handling thin structures, with or without finite losses (see Fig. 2) 


e Dealing with multi-scale geometries (see Fig. 3) 





Figure 1: An elliptical geometry. 





Figure 2: A very thin sheet. 


Advantages Some of the notable features of yFDTD are: 


e Unlike the conventional FDTD, the mesh-size utilized by the yFDTD is not dictated 


by the finest feature of the geometry, and this size is usually maintained at the 





Figure 3: PEC loop over a finite ground plane. 


conventional \/20 level, where \ corresponds to the highest freqeuncy of interest. 


This helps to reduce the computational burden by a large factor. 


e The vFDTD algorithm incorporates a novel post-processing technique which requires 
relatively few time steps, in comparison to the number of steps required by the 


conventional FDTD. 


0.3 Low Frequency Response 


Despite many advances in finite methods, such as the FEM (Finite Element Method) and 
the FDTD, as well as in integral-equation-based techniques such as the MoM (Method of 
Moments), it still remains a challenge to accurately calculate the low frequency response 
for radiation and scattering problems. The frequency domain techniques, such as the FEM 
and MoM, both experience difficulties at low frequencies, because they have to deal with 
ill-conditioned matrices at these frequencies. On the other hand, while the time-domain- 
based techniques, such as the FDTD, can accurately generate results at high frequencies, 


usually above 1 GHz, the same cannot be said about their performance at low frequencies. 


This is not only because the FDTD results are often corrupted by the presence of non- 
physical artifacts at low frequencies, but also because the FDTD requires an exorbitantly 
large number of time steps for accurate calculation of the response. The required number 
of time steps can exceed a few million in some cases before convergence is achieved. 

As an example, let us consider a 32 port connector circuit example shown in Fig.4. 
This connector geometry has been analyzed by using a commercial FDTD solver and the 
variation of the transmission coefficient S21 is plotted in Fig.5 as a function of the frequency, 
and we observe that the results show ripples that are numerical artifacts. Table 1 compares 
the number of time steps required for the solution to converge at different frequencies for 
the connector geometry. It can be inferred from this Table that the number of time 
steps required for the convergence increases as we go down in frequency, and eventually 
it becomes totally impractical to solve the problem at very low frequencies. Accurate 
calculation of the low frequency response becomes especially critical in the area of RF 
and digital circuits, since inaccurate results can affect the causality behavior of the overall 
system. The vFDTD utilizes a new technique, which is based on analytic continuation 
of the results derived at higher frequencies, and which is implemented by using the DM 
Approach and related techniques. This new technique is universal in nature, and it covers 
the entire range of frequencies, including the limiting case of f > 0. Also, the yFDTD can 
handle both the RF/Digital circuit problems as well as the radiation/scattering problems 
with same ease, by employing unique methodologies tailored for each of these categories. 


We present these methodologies in detail in the sections that follow. 


Table 1: Comparison of time steps required for convergence for the circuit shown in Fig.4. 





Frequency 10 MHz 1 MHz 1 Hz 
Time Steps in Millions | 0.7 7 70 








Figure 4: A 32 port connector with a overall dimension of 5.6 x 11.88 x 27.35 mm (housing 
not shown here). 
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Figure 5: Variation of the transmission coefficient S21 for the 32 port connector shown in 
Fig.4. 


0.3.1 RF and Digital Circuits 


Consider the variation of the isolation coefficient $3; shown in Fig.6 for the connector 


geometry(Fig.4). This plot is divided into three regions, namely: 


e Region-1: Low-frequency regime 
e Region-2: Validation region 


e Region-3: High-frequency regime 
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Figure 6: Variation of the isolation coefficient S31 for the 32 port connector shown in Fig.4. 


There are four frequency values which delimit the above three regions. The frequency 
fi describes the lowest frequency of interest defined by the user. The frequency fı, which 
divides the regions 1 and 2, is typically chosen to be between 500 MHz to 1000 MHz, while 
the frequency fə dividing the regions 2 and 3 is chosen to be on the order of 2f; or 3/1. 
The frequency fy is the user input indicating the highest frequency of interest. In each 
of these three regions the results are calculated by using a different method. The results 


in the high frequency regime are generated by using the conventional FDTD, using a DC 


Gaussian pulse as the excitation source, whose 3dB cut-off frequency is set to be fy. In 
the low frequency regime, the results are generated by using the proposed new technique, 


which involves the following steps: 


1. Smooth the “DC Gaussian” results. 


2. Fit the curve from fz to fı with the DC values, using a quadratic, for instance. The 


choice of fı can be fine-tuned based on the quality of the resulting fit. 


3. Validate the smoothed “DC Gaussian” results in region-2 by comparing them with 


those generated by “single frequency” simulations at a few points (typically 2 or 3). 


We have recalculated the results for the 32-port connector geometry, shown in Fig.4, by 
using the above method. The new results for the variation of the transmission coefficient 
So, and the isolation coefficient S31 are shown in Figs. 7 and 8. From these figures we 
can clearly see that the conventional FDTD simulation utilizing the DC Gaussian pulse 
as excitation does not generate an accurate low frequency response and has numerical 
artifacts, while the yFDTD does not suffer from the same. 

For the next example, we consider an 8-port connector as shown in Fig. 9, which 
operates in the frequency range 50-800 MHz. Fig. 10 compares the variation of S31 
calculated by using the vFDTD with those obtained by using the DC Gaussian in the 
conventional FDTD algorithm. Again we find that the DC Gaussian results show spurious 


spikes, while the yF DTD was able to calculate it accurately. 


0.3.2 Scattering Problems 


In this section we turn to the solution of scattering problems by using the vFDTD. The 
methodology for handling the radiation and scattering problems are different from those 


used for the RF/Digital circuits, as we will explain below. For the high frequency regime, 
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Figure 7: Variation of the transmission coefficient S21 for the 32 port connector shown in 
Fig.4 calculated using vFDTD. 





Magnitude in dB 
k 











sol = = = DC Gaussian 
Æ Single Frequency 
vFDTD 
-55 L 
-60 fi 1 J 
500 1000 1500 


Frequency in MHz 


Figure 8: Variation of the isolation coefficient S31 for the 32 port connector shown in Fig.4 
calculated using vFDTD. 
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Figure 10: Variation of the $3; for the 8 port connector shown in Fig.9. 


we use the conventional FDTD, and use a Gaussian excitation source to generate the 
results. However, we employ a different procedure, as outlined below, in the low frequency 


regime: 


1. Run a “Single Frequency” simulation at a frequency fı, where the largest dimension 


of the geometry is A/100, to calculate the fields at a point located 4/20 from the 
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surface of the object. 


2. Extract the dipole moment by using the analytical expressions for the field radiated 


by an infinitesimal dipole [4] from the field values calculated above. 


3. Use the extracted dipole moment to calculate the results from fzr to fo, where fzr is 
the lowest frequency of interest, and fə is typically chosen to be 2fı or 3fı. It has 
been found that the results generated by using this dipole moment is not only valid 
for frequencies as low as 0, but also up to frequencies where the largest dimension of 
the geometry becomes /10; hence it enables us to dovetail the low frequency results, 


seamlessly, with the lower end of the high frequency regime. 


4. Validate the “DC Gaussian” results in the range between fı and f2 by comparing 
them with those calculated by using the analytical expression at a few points (typi- 


cally 2 or 3). 


In order to extract the dipole moment from the single frequency simulation, one can 
either use the method proposed by Furse [5], or use the DFT to process the time signature. 
In the Furse method, we choose two samples of the time signature and we fit the time 
signature to a sinusoidal curve using those two samples. Even though this method looks 
computationally inexpensive when compared to the DFT approach, the choice of the two 
samples determines the accuracy of the method, and these samples should not lie within the 
transient region; hence we always use the DFT to extract the DM because of its robustness. 

As an example application of the procedure just outlined, we consider a sphere with a 
diameter of A/20, with À defined at 10 GHz. The sphere is illuminated by a plane wave 
traveling in the negative-z direction, with its E-field polarized along y. Fig.12 compares 
the fields calculated by the proposed technique, in the frequency range of 1 Hz to 30 GHz, 


with those derived analytically. We find that the fields calculated by using the proposed 
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technique based on DM extraction exhibits good agreement with those calculated by using 
the analytical expression. The small deviation between the two curves is attributable to 
the staircase modeling of the sphere in the conventional FDTD, and it can be corrected 
by using an effective radius in the analytical expression. It is important to recognize 
the fact that we have used the same technique to calculate the response over the entire 
frequency range, including frequencies as low as 1 Hz, without using either the quasi- 
static approximation or other special treatments that are employed in the conventional 
computational electromagnetic (CEM) techniques. Even after the use of special treatments 
in the existing techniques, such as the FEM and MoM, the accuracy of the low-frequency 
solution is often questionable because of the large condition number of the associated 
matrix. Thus, despite all the special treatments implemented in these methods to address 
the low frequency breakdown problem, it is totally impractical to go down to frequencies 
on the order of 1 Hz in the existing techniques. 

The amplitude variation of the scattered field with the distance along z, calculated by 
using the proposed technique, is shown in Fig.13 for a frequency of 1.8 GHz. This plot also 
compares the results with those calculated by using analytical expressions. Again we find 
good agreement between the yFDTD results and those generated from the analytical ex- 
pression for a \/67 sphere, for the chosen frequency of 1.8 GHz. The field variation derived 
by using the yFDTD matches well with that generated from the analytical expression, both 
in the near and far field regions. 

For the next example we consider a PEC cube of side length 1/20, as shown in Fig. 14 
at a frequency of 10 GHz. Fig. 15 plots the scattered E,-field as a function of frequency, 
calculated at z = 2.5mm (A/12 at 10 GHz), by using the vFDTD, and compares them 
with those obtained by using the DM approach. The comparison is seen to be good even 


at a frequency of 1 Hz. 
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Figure 11: A PEC sphere of diameter pa at 10 GHz. 
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Figure 12: Amplitude variation of the scattered E, at a point z = 0.25cm with frequencies 
from 1 Hz to 30 GHz. 


Based on the illustrative examples presented above, we list below some of the advantages 


of the proposed method: 
e RF and Digital Circuit Problems: 
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Figure 13: Amplitude variation of the scattered E, with distance along z from È to x, at 


1.8 GHz. 






=> Obs. Point : A/12 





A @ 10 GHz 


Figure 14: A PEC cube of side length A at 10 GHz. 


Efficient for constructing an accurate low frequency solution, compared to the long 


runs in FDTD. 


e Scattering Problems: 


(a) Can be used for an arbitrary geometry. 
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Figure 15: Amplitude variation of the scattered Ey at a point with frequencies from 1Hz 
to 30 GHz. 


(b) Can be used to efficiently calculate not only the frequency response, but the 


near and far fields as well. 


0.4 Non-Cartesian Geometries 


The conventional FDTD uses a staircase approximation to model non-Cartesian geometries, 
as shown in Fig.16, and requires the use of a very fine mesh to mitigate the effects of 
this staircase approximation when dealing with curved objects. This, in turn, makes the 
simulation computationally expensive, both in terms of memory and CPU time. Even 
though methods such as FEM and MoM can handle curved geometries with ease because 
they do not restrict themselves to a Cartesian type of meshing, often they are not necessarily 
the most computationally efficient when dealing with inhomogeneous media. Hence, it 


would be advantageous to modify the existing FDTD algorithm so that it can handle 
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curved geometries, enabling us to conveniently model arbitrary objects, regardless of their 
material parameters. In the past, a generalization of the conventional FDTD, namely the 
CFDTD algorithm [1], has been developed for this purpose. In CFDTD, the magnetic field 
update equations are modified by using the areas of the partially-filled cells, as opposed to 
those of the entire cells. 

To explain the concept, we consider a partially filled cell, shown in Fig. 17. The 


equation for this partially-filled cell is derived by using Farady’s law, to get: 











(b) A PEC wedge with staircase approx- 
(a) A PEC wedge geometry imation 


Figure 16: Meshing of a non-Cartesian geometry by the conventional FDTD 


-> > 0 => => 
E -dl = -p— H -ds 1 
f, mJ, a) 


where C is the loop ABCDA and S; is the area enclosed by the loop C1. Upon discretizing 


this equation, we obtain: 
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Figure 17: A partially-filled cell. 


m k) = ee k) — a Elid k) -lag + ER(i,j,k)-dh+ E? (i+1, j,k) lop] 
(2) 

The magnetic field update equation for the partially-filled cell is shown above in (2). 
But, as S; — 0, this modified update equation becomes unstable since, as we see from (2), 
the expression for the updated H contains S; in the denominator. The update equation 
can be modified to circumvent this instability problem that arises when the partial area 


is small, albeit at the cost of compromising the accuracy. Hence, in order to improve the 


accuracy, we propose two new approaches to handling the non-Cartesian geometries. 


0.4.1 Asymptotic Method 


In this asymptotic type of implementation, the field values, as opposed to the update 
equations, are modified by using the local field solution. The proposed new technique is 


described below: 


e For the partially-filled cells with a fill factor < 50%, the E-fields are updated by using 


the H-fields calculated by using the modified CFDTD equation given in (2). 
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e For the partially-filled cells with a fill factor > 50%, the E-fields are updated by using 
local solutions generated based on the concepts of reflection or diffraction, rather than 


by using the H-fields employed in the CFDTD approach. 


Because we use the asymptotic method to compute the reflection or diffraction coef- 
ficients, the proposed technique requires a “single frequency” simulation. However, this 
technique can be extended to “DC Gaussian” simulations with a slight modification as 
described in the next Section 0.4.2. Also, the proposed technique can be extended to 
dielectrics and inhomogeneous geometries without any modification, while the CFDTD 
cannot handle either of them without compromising the accuracy. 

Let us consider the case of a square PEC sheet whose sides are approximately 4\(A 
referenced at 10 GHz), and are inclined at an angle of 0.72° with respect to the x-axis, 


as shown in Fig.18. The tilt angle has been chosen to be 0.72° so that the edges of the 





sheet are offset only by +A/40 above or below the x-axis; i.e., half the FDTD cell size 
of A/20. We calculate the amplitude variation of the scattered EF, field at a frequency of 
10 GHz, when the plate is illuminated by a plane wave, which travels along the negative- 
y direction and is polarized along x. Fig.19 compares the results obtained by using the 
proposed technique, with those returned by the CFDTD and by a commercial MoM code, 
for the same problem. The results generated by using the proposed technique show good 
agreement with the results from the commercial MoM, while the CFDTD results exhibit 


spurious ripples in the lit region because of the instability problem encountered in the 


-À 
160 


CFDTD algorithm when the area S$; — 0. This is true even when a fine mesh size of is 
used in the CFDTD, in contrast to the pa mesh size used in vFDTD. Table 2, presents a 
comparison of the mesh size and the memory requirements, and shows that the proposed 
technique easily out-performs the CFDTD, which still suffers from inaccuracies, even when 


a very fine mesh is used. 
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Figure 18: A inclined PEC sheet(not to scale). 
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Figure 19: Amplitude variation of the scattered FE, with distance along y at 10 GHz. 
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Table 2: Comparison of mesh size and memory required for convergence for PEC geometry 
shown in Fig.18. 





Parameter | vFDTD CFDTD* 
Mesh Size Used a in 
Memory Required 413 MB 31 GB 





a Results still have numerical artifacts 


For the next example, let us change the inclination of the PEC plate in the previous 
problem from 0.72° to 1.43°. Fig. 20 compares the scattered E,-field calculated by using 
the vVFDTD/Asymptotic method, with those obtained by using the CFDTD algorithm with 
a mesh size of A/160, and with the commercial MoM code results. We again find that the 
results from the vFDTD match well with those obtained by using the commercial MoM, 
while the results from the CFDTD show spurious ripples due to the instability problem 


described above. 
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Figure 20: Amplitude variation of the scattered EF, with distance along y at 10 GHz. 
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For the next example, we consider a faceted PEC surface (see Fig. 21) whose projected 
length is À at a frequency of 10 GHz. Fig. 22 compares the backscattered E-field calculated 
by using the vFDTD with those obtained from: (i) the CFDTD with a mesh size of \/20; 
(ii) a commercial MoM code; and (iii) a commercial FEM code. Again we find that 
the results calculated by using the vFDTD compares well with those obtained from the 
commercial MoM code, while the results from the commercial FEM code show numerical 


artifacts. 


A @ 10 GHz 





Figure 21: A faceted PEC surface (not to scale). 


For the next example, let us consider a PEC wedge of side length 4X, as shown in Fig. 
23. Fig. 24 plots the scattered field at a frequency of 10 GHz along the specular direction, 
obtained by using yvFDTD, and compares them with those obtained by using the CFDTD 


with a mesh size of 4/50; with a commercial MoM code; and, with a commercial FEM 
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Figure 22: Amplitude variation of the backscattered Æ, with distance along y at 10 GHz. 


code. We find a good comparison between the scattered fields calculated by using the 
VFDTD with those obtained from the commercial MoM code, while the results from the 
CFDTD and the commercial FEM codes show spurious ripples. 

For our next example, we consider a finite PEC cylinder, as shown in the Fig. 25, with 
a height of 21/20 at a frequency of 10 GHz, for the next example. In order to calculate 
the field in the asymptotic limit in the yYFDTD method, we use the fields scattered by an 


infinite PEC sheet multiplied by the factor f, defined in (3) 


ae (3) 


where a’ is the effective phase center and r is the distance of the observation point from the 
effective phase center. A wide range of numerical experiments have shown that this effective 


phase center for a PEC cylinder is always 0.5a, where a is the radius of the cylinder. 
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Fz y A @ 10 GHz 
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Figure 23: A PEC wedge. 
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Figure 24: Amplitude variation of the scattered Æ, with distance along x at 10 GHz. 
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E, A @ 10 GHz 





Observation Line 


Figure 25: A PEC cylinder. 


Fig. 26 compares the backscattered E,-fields calculated by using the yFDTD with those 
obtained from the CFDTD with a mesh size of \/20; with a commercial MoM code; and, 
with a commercial FEM code. We find that the fields calculated by using the yFDTD shows 
good comparison with those obtained by using the commercial MoM solver. The results 
calculated using the commercial FEM solver shows numerical artifacts and the CFDTD 
does not generate the correct scattered field on the surface of the cylinder; however, the 
VFDTD is able to solve this problem with good accuracy. 

While modeling dielectric objects using the CFDTD approach, the partially-filled cells 
are smeared with an average dielectric constant over the entire volume of the cell. The 
asymptotic method proposed here can be used as an alternative to model dielectric objects 


without any modifications, and with better accuracy than that of the CFDTD method. As 
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Figure 26: Amplitude variation of the backscattered Æ, with distance along x at 10 GHz. 


an example, let us consider the dielectric slab of thickness \/8 at a frequency of 10 GHz, 
with er = 4.2 as shown in Fig. 27. Fig. 28 compares the backscattered E-field calculated 
by using the vFDTD with those generated by using the CFDTD, and with the commercial 
FEM code. We find that the results generated by using the commercial FEM and CFDTD 
codes show spurious ripples, while the those from the yFDTD have a smooth behavior, 
which is realistic. 

In order to demonstrate the efficacy of the vF DTD method, we vary the thickness of the 
slab, and choose it to be 9/80, 10/80 and 11X/80 at a frequency of 10 GHz. We calculate 
the phase variation at y = /40 and plot it in Fig. 29, which shows the comparison of 
the phase variation against the thickness, generated by using the infinite slab analytical 
expression; the CFDTD; and a commercial FEM code. We find that the results calculated 
by using the yFDTD compares well with those obtained from the analytical expressions 


for the infinite slab, while the results generated by using other methods deviate from the 
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Figure 27: A dielectric slab (not to scale). 
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Figure 28: Amplitude variation of the backscattered Æ, with distance along y at 10 GHz. 
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analytical results. 
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Figure 29: Phase variation of the backscattered E, at y = /40. 


0.4.2 Averaging Technique 


As mentioned earlier in Section 0.4.1, the proposed vFDTD technique requires a “single 
frequency” simulation because we use the asymptotic limit to compute the reflection or 
diffraction coefficients. In this section we describe a modified approach in which we use 
the “DC Gaussian” simulation that enables us to generate results for a wide range of 
frequencies. In order to demonstrate the usefulness of this approach, we consider the 
example of a rectangular cylinder of height 4A at a frequency of 10 GHz, which has mitered 
corners as shown in Fig. 30. The problem is handled in three steps as shown in Fig. 31 and 
we use scattered field type of formulation in all of these steps. Even though it shows three 
steps, it is important to note the fact that we need only two simulations as the calculation 


of the scattered field in step 2 is trivial using the boundary condition. The simulation for 
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the first step does not call for any modifications to the FDTD algorithm. However, for 
the third step we need to modify the field value at the nodes of the partially filled cells by 
using the weighted average of the fields at these nodes, obtained from steps 1 and 2, based 
on the partially filled space in the actual geometry. 
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Figure 30: A PEC rectangular cylinder with mitered corner. 


Fig. 32 compares the backscattered E,-fields calculated by using the vFDTD, and 
with those generated by using the commercial MoM code and the comparison is seen to 
be good. For the next example, we consider the faceted PEC geometry shown in Fig. 
33, whose projected footprint is 4A at a frequency of 10 GHz. Fig. 34 compares the 
backscattered E,-field, calculated by using the vFDTD/Averaging technique, and with 
those obtained from a commercial MoM solver. From Fig. 34, we find that the solution 


generated by using the commercial MoM code shows a spurious spike near the surface of 


28 





E1 (1) 


E,(1) 











E',(2) 


E,(2) 








E',(3) 


E,(3) 








E',(4) 





E,(4) 








E'\(5) 








E,(5) 























Actual Geometry 
Step 3 (A/20) 


Step 1 (A/20 Simulation) Step 2 (0, simulation not required) 


Figure 31: Principle behind the averaging technique. 
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Figure 32: Amplitude variation of the backscattered E, with the distance along x at 10 
GHz. 
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the geometry, while the yFDTD results are smooth. 
A @ 10 GHz 


y E 





Figure 33: A faceted PEC geometry (not to scale). 


For the last example, we consider a curved PEC surface of height 4 at a frequency 
of 10 GHz, as shown in Fig. 35. Fig. 36 compares the scattered E,-fields calculated by 
using the yFDTD method, and with those from a commercial MoM code, and once again 


we find that the comparison is good. 


0.4.3 Advantages 


Below, we summarize some of the advantages of the proposed method. They are: 


(a) Usable for arbitrary geometries, even if the surfaces do not coincide with the Cartesian 


mesh, e.g., thin sheets, with or without a slant. 


(b) More accurate than the conventional conformal FDTD. 
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Figure 34: Amplitude variation of the backscattered Ey with distance along x at 10 GHz. 


A @10GHz y 


7N/20 





Observation Line 





<< 
A/10 


Figure 35: A curved PEC surface (with a height of 4A). 


(c) Retains \/20 cell size even for thin, slanted and curved bodies, offering memory ad- 


vantage and computational efficiency over conventional conformal FDTD. 
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Figure 36: Amplitude variation of the scattered Ey at 10 GHz. 


(d) Free of instability problems even when the fractional area of the partially-filled cell is 


very small, even when it tends to zero. 


(e) Can be extended to dielectric objects, with just a few modifications. 


0.5 Multiscale Problems 


The direct solution of multiscale problems by means of conventional CEM methods—be 
it FEM, FDTD or MoM-is highly challenging, even with the availability of modern su- 
percomputers, because we need to use a large number of DoFs (degrees of freedom) to 
accurately describe objects with fine features, which might share the computational do- 
main with other large objects. For instance, if a thin wire is located in an inhomogeneous 
medium along with other large-scale objects, and the thickness of the wire is only a small 


fraction of the wavelength in the medium, then we must use a very fine mesh to accurately 
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capture the nuances of its geometry. This, in turn, leads to a large number of DoFs when 
other large objects are also present in the computational domain along with the thin wire 
whose shape may be arbitrary. 

Dealing with multiscale objects often forces us to compromise the accuracy (relaxing 
the numerical discretization process when attempting to capture the small-scale features) 
in order to cope with the limited available resources in terms of CPU memory and time. 
In this section, we introduce a scheme that combines the FDTD and the DM approach 
to solve multiscale problems in a numerically efficient manner. Our objective is to handle 
objects with fine features with the DM approach, and not directly with the FDTD which 
would require us to use a fine mesh (see Fig. 38 for the problem shown in Fig. 37), at the 
cost of an increased computational burden when compared to that of a problem without 
fine features. 

The main advantage of this hybrid method is that it does not require local mesh refine- 
ment for objects with fine features (Fig. 39). In fact, the region surrounding the small/thin 
structure is extracted from the original domain and two different numerical techniques are 
used for dealing with the two problems. The coupling of the object with the remaining 
part of the computational domain is achieved by using the fields radiated by the previously 
extracted region. As a result, the presented method does not place a heavy burden on the 
CPU time and memory as do the conventional approaches when dealing with multiscale 
problems. The DM/FDTD method introduced herein can be implemented either in an 
iterative or in a self-consistent manner. 

The proposed method is especially useful for modeling wire antennas located in the 
vicinity of inhomogeneous structures, as well as for simulating interconnect structures in 


integrated circuits, which typically have fine features. 
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Figure 37: A multiscale problem. 


0.5.1 Iterative Approach 


Both the iterative and self-consistent hybrid implementations - the latter to be described 
in the following section - begin by extracting a region surrounding the small object from 
the FDTD domain. A 2-D representation of the hybrid problem is shown Fig. 39. 

Let us assume that two objects, a large PEC plate and a PEC wire, which is small 
compared to the operating wavelength, are located in the FDTD computational domain, 
which is excited by a voltage gap source, at 10 GHz as shown in Fig. 40. The hybrid- 
iteration algorithm begins by solving the small object (dipole antenna in this case), which 


may be PEC, or dielectric (or a combination thereof) and is treated by using the DM 
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With PEC Object 


Figure 38: A multiscale problem meshed for FDTD simulation. 


approach described in [6], in the absence of the large structure. 

Next, the fields scattered by the large structure are derived by using FDTD and a 
source excitation comprising of the fields radiated by the small object. 

These scattered fields are evaluated at the boundary of the extracted region and then 
interpolated to obtain the incident fields at the locations of each of the basis functions used 
to model the dipole in DM approach. Following this, the right hand side of the matrix 
equation, generated using the DM approach for the dipole is modified, by superimposing 


the feed and the fields scattered by the larger structure. Next, the matrix equation is 
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Figure 39: A multiscale problem meshed for hybrid FDTD simulation. 


solved for the weight coefficients of the dipole moments as a first step in the iteration 
process. Then, we again solve for the fields scattered by the large object when illuminated 
by the fields radiated by the small object, derived by using the recently calculated weight 
coefficients of the current in the small object. 

The iteration process is continued, by repeating the steps described above. The process 
is terminated when numerical convergence has been achieved and the difference between 
the results obtained at the k-th and (k-1)-th iteration steps is below a chosen threshold, 


say, 1073. Fig. 41 compares the amplitude of the scattered field Ey calculated by using 
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Figure 40: A A/20 dipole antenna over a finite ground plane (not to scale). 


the FDTD iterative hybrid approach with those obtained from a commercial MoM solver. 

Fig. 41 shows that the scattered fields calculated by using the hybrid FDTD iterative 
method compares well with those obtained from a commercial MoM code. The commercial 
MoM code was able to handle this problem with ease since the large object was modeled 
as a PEC sheet; however, when the PEC sheet is replaced by a thin PEC plate, the CPU 
time and the memory required by the commercial MoM increases, as will be demonstrated 


in the next section. 
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Figure 41: Amplitude comparison of Ey field for the multiscale problem shown in Fig. 40. 


0.5.2 Self-Consistent Approach 


The self-consistent type of hybrid implementation also begins by extracting from the FDTD 
domain a region surrounding the small object (see Fig. 39). However, this time, the entire 
problem is solved in a single step by directly inverting a composite matrix equation, which 
is constructed as follows: First, the impedance matrix for the small problem is set up 
independently of the rest, by using the DM approach; the right-hand side vector for source 
is computed and stored. Next, we compute the field radiated by the current distribution 
on the small object at the location of the large object, and solve for the scattered field 
on its surface by imposing the boundary condition with the fields produced by the small 


object as the incident field on the large object. The fields scattered by the large object 
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are computed in the entire computational domain by using a FDTD simulation carried out 
by using a coarse mesh, and then interpolated in the region containing the small object 
to obtain a new excitation vector for the DM system. Now the matrix equation is solved 
for the weight coefficients associated with the currents on the small object, and the final 
scattered fields are calculated as a weighted superposition of the contributions from the 
small and large objects. 

To illustrate the procedure, let us consider a dipole antenna of length /2 over a finite 
ground plane, operating at 10 GHz as shown in Fig. 42. Fig. 43 compares the scattered field 
calculated by using the self-consistent approach as described above, with those calculated 
by using the iterative approach and the commercial MoM code. Table 3 compares the 
computational resources required by the self-consistent and iterative approaches. 


Table 3: Comparison of computational resources required by the different FDTD hybrid 
approaches for the multiscale problem shown in Fig.42. 





| Self-Consistent Approach Iterative Approach 
Peak Memory 483 MB 481 MB 
Simulation Time | 172.3 s 601 s 





Fig. 43 shows good comparison of the scattered field calculated by using the self consis- 
tent approach, with those calculated by using the iterative procedure and the commercial 
MoM code, with the exception of observation points located near the source region. Ta- 
ble 3 shows that the results calculated by using the iterative approach are more accurate, 
though they have a longer run time in comparison to the self-consistent approach. In order 
to demonstrate the key advantage of using the hybrid FDTD, the PEC sheet in Fig. 42 
was replaced by a PEC plate of thickness 4/20. The problem was solved by using the 
self-consistent approach and the variation of scattered field E, is plotted in Fig. 44. Ta- 


ble 4 compares the computational resources required by the self-consistent approach with 


39 






Feed Gap Source E, (A/20) 
1.5A/20 
13A/20 


Figure 42: A \/2 dipole antenna over a finite ground plane (not to scale). 


those needed by the commercial MoM. From Fig. 44 and the Table 4, we find that the 
self-consistent approach outperforms the commercial MoM code, both in terms of memory 
and simulation time. 


Table 4: Comparison of computational resources required by the self-consistent approach 
with those required by the commercial MoM. 





oe Hy- Comm. MoM 


brid 
Peak Memory 248 MB 543 MB 
Simulation Time 71s 224.75 s 
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Figure 43: Amplitude comparison of Ey field for the multiscale problem shown in Fig. 42. 


0.6 Observations and Conclusions 


In this write-up, we have introduced the yFDTD solver, which is a blend of time and 
frequency domain techniques that can generate accurate electromagnetic responses at low 
frequencies; handle non-Cartesian geometries accurately without any instability issues that 
are often encountered in the conventional CFDTD; model multi-scale geometries accu- 
rately; and handle lossy/lossless thin structures with ease. In all the cases for which we 
have carried out comparison studies with the existing algorithms and commercial codes, 
the vFDTD was not only accurate but also computationally the most efficient. Finally, 
since yvFDTD builds on the conventional FDTD to solve different types of problems, its per- 


formance can be further enhanced by parallelizing the algorithm [7], which can be carried 
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Figure 44: Amplitude comparison of Ey field for the multiscale problem with a \/20 thick 
finite ground plane. 


out easily as in the case of the conventional FDTD. 
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